function [A, b] = assemb_system(dof_map, K11, K12, K22, b1, b2, fb_1, fb_2)
%function [A, b] = assemb_system(dof_map, K11, K12, K22, b1, b2, fb_1, fb_2)
%
%

% firstly, assemble element system into large sparse matrix
n_elem = size(K11, 3);
n_dof_per_elem = size(K11, 1);
n_nz = n_dof_per_elem*n_dof_per_elem;
n_dof = max(max(dof_map));
% n_bezier = size(dof_map,1)*size(dof_map,2);

i = zeros(n_nz, n_elem);
j = i;
K11s = reshape(K11, n_nz, n_elem);
K12s = reshape(K12, n_nz, n_elem);
K22s = reshape(K22, n_nz, n_elem);

% b = zeros(2*n_dof, 1);
b = [fb_1; fb_2]; % initially the surface force

i_tmp = repmat((1:n_dof_per_elem)', n_dof_per_elem, 1);
j_tmp = reshape(repmat((1:n_dof_per_elem) , n_dof_per_elem, 1), n_nz, 1);
for tri = 1:n_elem
     i(:, tri) = dof_map(i_tmp, tri);
     j(:, tri) = dof_map(j_tmp, tri);
     b(dof_map(:,tri)) = b(dof_map(:,tri)) + b1(:, tri);
     b(dof_map(:,tri)+n_dof) = b(dof_map(:,tri)+n_dof) + b2(:, tri);
end

A = sparse(i, j, K11s, 2*n_dof, 2*n_dof) + sparse(i, j+n_dof, K12s, 2*n_dof, 2*n_dof) + ...
       sparse(j+n_dof, i, K12s, 2*n_dof, 2*n_dof) + sparse(i+n_dof, j+n_dof, K22s, 2*n_dof, 2*n_dof);

   
% b = [reshape(b1, n_bezier, 1) ; reshape(b2, n_bezier, 1)];

% secondly,  applied the dirichlet boundary
% if (nargin > 6)
%     [n_row, n_col] = size(dof_boundary);
%     for j = 1: n_col
%         for i = 1: n_row
%             dof = dof(boundary(i, j)); 
%             A(dof,:) = 0; A(dof, dof) = 1;
%             A(n_dof + dof,:) = 0; A(n_dof + dof, n_dof + dof) = 1;
%             b(dof) = ub_1(i,j); 
%             b(n_dof + dof) = ub_2(i,j);
%             b = 
%         end
%     end
% end

end